Exploración de datos

Modificamos los datos para ver las medias por artista:

spotify <- spotify %>% 
  select(artist, title, popularity) %>% 
  mutate(artist = fct_reorder(artist, popularity, .fun = 'mean'))

artist_means <- spotify %>% 
  group_by(artist) %>% 
  summarize(count = n(), popularity = mean(popularity))

Los artistas más y menos populares:

artist_means %>%
  slice(1:2, 43:44)
## # A tibble: 4 × 3
##   artist        count popularity
##   <fct>         <int>      <dbl>
## 1 Mia X             4       13.2
## 2 Chris Goldarg    10       16.4
## 3 Lil Skies         3       79.3
## 4 Camilo            9       81

Modelo full pooleado

Tenemos de 2 a 40 canciones por artista:

artist_means %>% 
  summarize(min(count), max(count))
## # A tibble: 1 × 2
##   `min(count)` `max(count)`
##          <int>        <int>
## 1            2           40

La estimación no paramétrica de la densidad tiene esta pinta:

ggplot(spotify, aes(x = popularity)) + 
  geom_density(linewidth = 1, fill = "steelblue", alpha = .3) +
  theme_bw()

Ahora vamos a ajustar el modelo full pooleado:

spotify_complete_pooled <- stan_glm(
  popularity ~ 1, 
  data = spotify, family = gaussian, 
  prior_intercept = normal(50, 2.5, autoscale = TRUE),
  prior_aux = exponential(1, autoscale = TRUE),
  chains = 4, iter = 5000*2, seed = 84735)

Si bien en el texto dice que el prior de \(\mu\) es \(\mu \sim N(50, 52^2)\), podemos ver que en el código dice normal(50, 2.5, autoscale = TRUE). Esto es porque el autor aprovecha la opción autoscale de stan_glm(). En el caso del prior gaussiano esto significa que \(\sigma\) es \(2.5 \times S_y\):

round(2.5 * sd(spotify$popularity))
## [1] 52

prior_aux por defecto si el prior de intercept es normal es \(\sigma\). Para setear el prior_aux también aprovecha la opción autoscale de stan_glm(). Esto significa, para el caso exponencial, que el prior para \(\sigma\) es \(1 / S_y\):

round(1/sd(spotify$popularity), digits = 3)
## [1] 0.048

De hecho, si chequeamos los priors podemos ver los valore ajustados:

# Get prior specifications
prior_summary(spotify_complete_pooled)
## Priors for model 'spotify_complete_pooled' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 50, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 50, scale = 52)
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.048)
## ------
## See help('prior_summary.stanreg') for more details

Y ver el posterior summary:

complete_summary <- tidy(spotify_complete_pooled, 
                         effects = c("fixed", "aux"), 
                         conf.int = TRUE, conf.level = 0.80)
complete_summary
## # A tibble: 3 × 5
##   term        estimate std.error conf.low conf.high
##   <chr>          <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)     58.4     1.10      57.0      59.8
## 2 sigma           20.7     0.776     19.7      21.7
## 3 mean_PPD        58.4     1.57      56.4      60.4

Las estimaciónes del modelo en relación a las medias por artista:

set.seed(84735)
predictions_complete <- posterior_predict(spotify_complete_pooled,
                                          newdata = artist_means)

ppc_intervals(artist_means$popularity, yrep = predictions_complete,
              prob_outer = 0.80) +
  ggplot2::scale_x_continuous(labels = artist_means$artist,
                              breaks = 1:nrow(artist_means)) +
  theme_bw() +
  xaxis_text(angle = 90, hjust = 1)  +
  labs(x = "Artista", y = "Popularidad")

Las estimaciones promedio por artista van a ser todas iguales e iguales a las estimaciones por canción, ya que todas las estimaciones son iguales.

\[ E(\mu|y) \approx 58.39 \]

Modelo sin poolear

Ahora vamos a conservar todos los artistas por separado.

ggplot(spotify, aes(x = popularity, group = artist)) + 
  geom_density(linewidth = .5, fill = "steelblue", alpha = .3) +
  theme_bw()

\[ Y_{ij}|\mu_j,\sigma \sim N(\mu, \sigma^2) \]

Ahora definimos el modelo. Lo único que cambia es que agregamos artist como predictor (y con el -1 le sacamos el intercept). Fijensé que ahora en lugar de prior_intercept dice sólo prior y por eso es que en la definici’on del libro dice que \(\mu_j \sim N(50, S^2_j)\).

spotify_no_pooled <- stan_glm(
  popularity ~ artist - 1, 
  data = spotify, family = gaussian, 
  prior = normal(50, 2.5, autoscale = TRUE),
  prior_aux = exponential(1, autoscale = TRUE),
  chains = 4, iter = 5000*2, seed = 84735)

Acá también podemos chequear los priors:

# Get prior specifications
prior_summary(spotify_no_pooled)
## Priors for model 'spotify_no_pooled' 
## ------
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [50,50,50,...], scale = [2.5,2.5,2.5,...])
##   Adjusted prior:
##     ~ normal(location = [50,50,50,...], scale = [484.78,309.30,434.23,...])
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.048)
## ------
## See help('prior_summary.stanreg') for more details

Y ahora podemos ver los ajustes por artista con respecto a la \(\bar{y}_j\):

# Simulate the posterior predictive models
set.seed(84735)
predictions_no <- posterior_predict(
  spotify_no_pooled, newdata = artist_means)
  
# Plot the posterior predictive intervals
ppc_intervals(artist_means$popularity, yrep = predictions_no, 
              prob_outer = 0.80) +
  ggplot2::scale_x_continuous(labels = artist_means$artist, 
                              breaks = 1:nrow(artist_means)) +
  theme_bw() +
  geom_point(x = 1:nrow(artist_means), y = complete_summary$estimate[1], color = "darkorange") +
  xaxis_text(angle = 90, hjust = 1) +
  labs(x = "Artista", y = "Popularidad")

Qué pasa si poneoms el mismo prior para \(\mu\) pero con menor dispersión. Pongamos \(0.1\) en la función que sería equivalente a decir que \(\mu_j \sim N(50, 2)\):

spotify_no_pooled_smallsigma <- stan_glm(
  popularity ~ artist - 1, 
  data = spotify, family = gaussian, 
  prior = normal(50, .1, autoscale = TRUE),
  prior_aux = exponential(1, autoscale = TRUE),
  chains = 4, iter = 5000*2, seed = 84735)
# Simulate the posterior predictive models
set.seed(84735)
predictions_no <- posterior_predict(
  spotify_no_pooled_smallsigma, newdata = artist_means)
  
# Plot the posterior predictive intervals
ppc_intervals(artist_means$popularity, yrep = predictions_no, 
              prob_outer = 0.80) +
  ggplot2::scale_x_continuous(labels = artist_means$artist, 
                              breaks = 1:nrow(artist_means)) +
  theme_bw() +
  geom_point(x = 1:nrow(artist_means), y = complete_summary$estimate[1], color = "darkorange") +
  xaxis_text(angle = 90, hjust = 1) +
  labs(x = "Artista", y = "Popularidad")

Efectivamente vemos que los \(\mu_j\) quedan más cerca de \(50\) y ya no son la media de las observaciones para ese artista.

El modelo jerárquico

Ahora tenemos básicamente tres capas de modelado:

\[ \begin{array} _Capa 1: Y_{ij}|\mu_j,\sigma_y &\sim& Modela \,\, cómo \,\, varía \,\, la \,\, popularidad \,\, dentro \,\, del \,\, artista_j \\ Capa 2: \mu_j|\mu,\sigma_\mu &\sim& Modela \,\, cómo \,\, varía \,\, la \,\, popularidad \,\, \mu_j \,\, entre \,\, artistas \\ Capa 3: \mu,\sigma_y,\sigma_\mu &\sim& Priors \,\, para \,\, los \,\, parámetros \,\, globales \,\, compartidos \end{array} \]

En la definición formal del Christensen:

\[ \begin{array} _y|\theta,\alpha &\sim& f(y|\theta,\alpha) \\ \theta|\alpha &\sim& p_0(\theta|\alpha) \\ \alpha &\sim& p_1(\alpha) \equiv p_1(\alpha|\beta_0)\\ \end{array} \]

Podemos ver la densidad de las medias por artista:

ggplot(artist_means, aes(x = popularity)) + 
  geom_density(linewidth = 1, fill = "steelblue", alpha = .3) +
  theme_bw()

Ajustamos el modelos jerárquico. El término (1 | artist) indica que artista es un factor de agrupamiento (o efecto aleatorio). prior_covariance tiene que ver con la matriz de covarianza de los efectos aleatorios, sería el prior de \(\mu_j\) y, no entiendo bien por qué, sería lo mismo que poner un prior \(\mathcal{E}(1)\):

spotify_hierarchical <- stan_glmer(
  popularity ~ (1 | artist), 
  data = spotify, family = gaussian,
  prior_intercept = normal(50, 2.5, autoscale = TRUE),
  prior_aux = exponential(1, autoscale = TRUE),
  prior_covariance = decov(reg = 1, conc = 1, shape = 1, scale = 1),
  chains = 4, iter = 5000*2, seed = 84735)

Y vemos los priors:

# Confirm the prior tunings
prior_summary(spotify_hierarchical)
## Priors for model 'spotify_hierarchical' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 50, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 50, scale = 52)
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.048)
## 
## Covariance
##  ~ decov(reg. = 1, conc. = 1, shape = 1, scale = 1)
## ------
## See help('prior_summary.stanreg') for more details

Podemos ver el diagnóstico sobre los 47 parámetros (44 \(\mu_j\), \(\mu\), \(\sigma_y\) y \(\sigma_\mu\)):

mcmc_trace(spotify_hierarchical) + theme_bw()

mcmc_dens_overlay(spotify_hierarchical) + theme_bw()

neff_ratio(spotify_hierarchical)
##                              (Intercept) 
##                                  0.17855 
##              b[(Intercept) artist:Mia_X] 
##                                  0.77280 
##      b[(Intercept) artist:Chris_Goldarg] 
##                                  0.51375 
##          b[(Intercept) artist:Soul&Roll] 
##                                  0.66455 
##         b[(Intercept) artist:Honeywagon] 
##                                  0.79880 
##           b[(Intercept) artist:Röyksopp] 
##                                  0.75600 
##          b[(Intercept) artist:Freestyle] 
##                                  0.80445 
##           b[(Intercept) artist:DA_Image] 
##                                  0.81100 
##          b[(Intercept) artist:Jean_Juan] 
##                                  0.64785 
##           b[(Intercept) artist:TV_Noise] 
##                                  0.41910 
##          b[(Intercept) artist:Kid_Frost] 
##                                  0.83115 
##      b[(Intercept) artist:Tamar_Braxton] 
##                                  0.84595 
##              b[(Intercept) artist:BUNT.] 
##                                  0.87730 
##                b[(Intercept) artist:MTK] 
##                                  0.74510 
##       b[(Intercept) artist:Atlas_Genius] 
##                                  0.68040 
##           b[(Intercept) artist:Jazzinuf] 
##                                  0.62940 
##              b[(Intercept) artist:Elisa] 
##                                  0.79290 
##      b[(Intercept) artist:House_Of_Pain] 
##                                  0.68000 
## b[(Intercept) artist:Black_Stone_Cherry] 
##                                  0.65250 
##              b[(Intercept) artist:C-Kan] 
##                                  0.77330 
##          b[(Intercept) artist:Zeds_Dead] 
##                                  0.68750 
##     b[(Intercept) artist:David_Lee_Roth] 
##                                  0.86805 
##               b[(Intercept) artist:NODE] 
##                                  0.67540 
##   b[(Intercept) artist:Michael_Kiwanuka] 
##                                  0.80565 
##         b[(Intercept) artist:The_Wrecks] 
##                                  0.85610 
##             b[(Intercept) artist:The_xx] 
##                                  0.51225 
##            b[(Intercept) artist:Placebo] 
##                                  0.67910 
##  b[(Intercept) artist:Mike_WiLL_Made-It] 
##                                  0.63020 
##     b[(Intercept) artist:Kendrick_Lamar] 
##                                  0.31175 
##      b[(Intercept) artist:X_Ambassadors] 
##                                  0.57580 
##             b[(Intercept) artist:Hinder] 
##                                  0.62940 
##              b[(Intercept) artist:Au/Ra] 
##                                  0.66690 
##      b[(Intercept) artist:Missy_Elliott] 
##                                  0.70295 
##          b[(Intercept) artist:The_Blaze] 
##                                  0.81110 
##    b[(Intercept) artist:Vampire_Weekend] 
##                                  0.69715 
##               b[(Intercept) artist:Alok] 
##                                  0.37185 
##      b[(Intercept) artist:León_Larregui] 
##                                  0.83365 
##     b[(Intercept) artist:Sufjan_Stevens] 
##                                  0.81400 
##            b[(Intercept) artist:Beyoncé] 
##                                  0.34905 
##        b[(Intercept) artist:Frank_Ocean] 
##                                  0.28955 
##      b[(Intercept) artist:Sean_Kingston] 
##                                  0.93845 
##            b[(Intercept) artist:J._Cole] 
##                                  0.49660 
##     b[(Intercept) artist:Camila_Cabello] 
##                                  0.30635 
##          b[(Intercept) artist:Lil_Skies] 
##                                  0.84850 
##             b[(Intercept) artist:Camilo] 
##                                  0.58210 
##                                    sigma 
##                                  0.84050 
##    Sigma[artist:(Intercept),(Intercept)] 
##                                  0.19085
rhat(spotify_hierarchical)
##                              (Intercept) 
##                                1.0008220 
##              b[(Intercept) artist:Mia_X] 
##                                1.0001391 
##      b[(Intercept) artist:Chris_Goldarg] 
##                                0.9998651 
##          b[(Intercept) artist:Soul&Roll] 
##                                1.0002406 
##         b[(Intercept) artist:Honeywagon] 
##                                1.0001168 
##           b[(Intercept) artist:Röyksopp] 
##                                1.0000496 
##          b[(Intercept) artist:Freestyle] 
##                                1.0000014 
##           b[(Intercept) artist:DA_Image] 
##                                1.0000454 
##          b[(Intercept) artist:Jean_Juan] 
##                                1.0000395 
##           b[(Intercept) artist:TV_Noise] 
##                                1.0000438 
##          b[(Intercept) artist:Kid_Frost] 
##                                1.0000608 
##      b[(Intercept) artist:Tamar_Braxton] 
##                                0.9999553 
##              b[(Intercept) artist:BUNT.] 
##                                0.9999160 
##                b[(Intercept) artist:MTK] 
##                                1.0005003 
##       b[(Intercept) artist:Atlas_Genius] 
##                                1.0001059 
##           b[(Intercept) artist:Jazzinuf] 
##                                0.9999770 
##              b[(Intercept) artist:Elisa] 
##                                0.9999609 
##      b[(Intercept) artist:House_Of_Pain] 
##                                1.0000742 
## b[(Intercept) artist:Black_Stone_Cherry] 
##                                0.9999885 
##              b[(Intercept) artist:C-Kan] 
##                                0.9999436 
##          b[(Intercept) artist:Zeds_Dead] 
##                                0.9999352 
##     b[(Intercept) artist:David_Lee_Roth] 
##                                1.0001916 
##               b[(Intercept) artist:NODE] 
##                                1.0000100 
##   b[(Intercept) artist:Michael_Kiwanuka] 
##                                1.0000025 
##         b[(Intercept) artist:The_Wrecks] 
##                                0.9999947 
##             b[(Intercept) artist:The_xx] 
##                                1.0000460 
##            b[(Intercept) artist:Placebo] 
##                                1.0000410 
##  b[(Intercept) artist:Mike_WiLL_Made-It] 
##                                1.0000320 
##     b[(Intercept) artist:Kendrick_Lamar] 
##                                1.0004022 
##      b[(Intercept) artist:X_Ambassadors] 
##                                1.0001418 
##             b[(Intercept) artist:Hinder] 
##                                1.0002135 
##              b[(Intercept) artist:Au/Ra] 
##                                1.0000875 
##      b[(Intercept) artist:Missy_Elliott] 
##                                1.0000266 
##          b[(Intercept) artist:The_Blaze] 
##                                0.9999593 
##    b[(Intercept) artist:Vampire_Weekend] 
##                                0.9998935 
##               b[(Intercept) artist:Alok] 
##                                1.0003676 
##      b[(Intercept) artist:León_Larregui] 
##                                0.9999028 
##     b[(Intercept) artist:Sufjan_Stevens] 
##                                1.0000550 
##            b[(Intercept) artist:Beyoncé] 
##                                1.0002370 
##        b[(Intercept) artist:Frank_Ocean] 
##                                1.0006008 
##      b[(Intercept) artist:Sean_Kingston] 
##                                1.0000611 
##            b[(Intercept) artist:J._Cole] 
##                                1.0001423 
##     b[(Intercept) artist:Camila_Cabello] 
##                                1.0005154 
##          b[(Intercept) artist:Lil_Skies] 
##                                0.9999641 
##             b[(Intercept) artist:Camilo] 
##                                1.0000691 
##                                    sigma 
##                                1.0001962 
##    Sigma[artist:(Intercept),(Intercept)] 
##                                1.0007369

Y también la posterior:

pp_check(spotify_hierarchical) + 
  xlab("Popularidad") +
  theme_bw()

Podemos ver también las posteriores por artista:

set.seed(84735)
predictions_hierarchical <- posterior_predict(spotify_hierarchical, 
                                              newdata = artist_means)

# Posterior predictive plots
ppc_intervals(artist_means$popularity, yrep = predictions_hierarchical, 
              prob_outer = 0.80) +
  ggplot2::scale_x_continuous(labels = artist_means$artist, 
                              breaks = 1:nrow(artist_means)) +
  theme_bw() +
  xaxis_text(angle = 90, hjust = 1) + 
  geom_hline(yintercept = 52.5, color = "darkorange")  +
  labs(x = "Artista", y = "Popularidad")

artist_summary <- tidy(spotify_hierarchical, effects = "ran_vals", 
                       conf.int = TRUE, conf.level = 0.80)

Está interesante ver el tema de las variabilidades.

\[ \begin{array} _\frac{\sigma_y^2}{\sigma_\mu^2+\sigma_y^2} \\ \frac{\sigma_\mu^2}{\sigma_\mu^2+\sigma_y^2} \end{array} \]

artist_summary <- tidy(spotify_hierarchical, effects = "ran_vals", 
                       conf.int = TRUE, conf.level = 0.80)

artist_chains <- spotify_hierarchical %>%
  spread_draws(`(Intercept)`, b[,artist]) %>% 
  mutate(mu_j = `(Intercept)` + b) 

artist_summary_scaled <- artist_chains %>% 
  select(-`(Intercept)`, -b) %>% 
  mean_qi(.width = 0.80) %>% 
  mutate(artist = fct_reorder(artist, mu_j))

ggplot(artist_summary_scaled, 
       aes(x = artist, y = mu_j, ymin = .lower, ymax = .upper)) +
  geom_pointrange() +
  theme_bw() +
  xaxis_text(angle = 90, hjust = 1) +
  labs(x = "Artista", y = "mu_j")

set.seed(84735)
prediction_shortcut <- posterior_predict(
  spotify_hierarchical,
  newdata = data.frame(artist = c("Frank Ocean", "Mohsen Beats")))

# Posterior predictive model plots
mcmc_areas(prediction_shortcut, prob = 0.8) +
  ggplot2::scale_y_discrete(labels = c("Frank Ocean", "Mohsen Beats")) + theme_bw() + labs(x =  "Popularidad")
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.